Manipulation of heat current by the interface between graphene and white graphene 
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We investigate the heat current flowing across the interface between graphene and hexagonal boron 
nitride (so-called white graphene) using both molecular dynamics simulation and nonequilibrium 
Green's function approaches. These two distinct methods discover the same phenomena that the heat 
current is reduced linearly with increasing number of carbon atom at the interface, and the zigzag 
interface causes stronger reduction of heat current than the armchair interface. These phenomena are 
interpreted by both the lattice dynamics analysis and the transmission function explanation, which 
both reveal that the localized phonon modes at interfaces are responsible for the heat management. 
The room temperature interface thermal resistance is about 7 x 10~^°m'^K/W in zigzag interface 
and 3.5 x 10~^''m^K/W in armchair interface, which directly results in stronger heat reduction in 
zigzag interface. Our theoretical results provide a specific route for experimentalists to control the 
heat transport in the graphene and hexagonal boron nitride compound through shaping the interface 
between these two materials. 

PACS numbers: 65.80.-g, 63.22.-m, 61.48.-c, 68.35. Ja, 79.60.Jv 
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I. INTRODUCTION 



The one-atom-thick graphene sheet possesses unique 
electronic properties based on the zero energy gap in its 
electron band structure. However, a finite band gap is 
of crucial importance for applications in field-effect tran- 
sistor devices. Therefore, a meaningful work is to real- 
ize band gap opening in the electron band structure of 
graphene. Both experimentali and theoretical2iS groups 
have recently demonstrated a novel approach to accom- 
plish tunable band gap through growing hexagonal boron 
nitride islands within graphene sheet. This improve- 
ment benefits the practical development of graphene- 
based electronic nano-devices. For high qulity electronic 
devices, the Joule heating induced breakdown is a sig- 
nificant bottleneck to the advancement of the devicesr'"— 
so it is of essential importance to enhance the capability 
of heat delivery in these devices. The situation becomes 
completely opposite in the thermoelectric field, where low 
thermal conductivity is in great demand for high value of 
ZT.^ Hence, it is an urgent and practical task to supply 
some information for experimentalists on the modifica- 
tion of the heat transport capability for their graphene 
and hexagonal boron nitride (CBN) compound for differ- 
ent application motivations. That is the aim of present 
work. 

In this work, we use both molecular dynamics (MD) 
simulation and nonequilibrium Green's function (NEGF) 
approaches to investigate the heat transport across the 
CBN interface. The interface is built inner a graphene 
nanoribbon of armchair edge. We study the heat current 
reduced by CBN interface of different shapes and sizes, 
including zigzag trigonal interface [CBN(TZ)], armchair 
trigonal interface [CBN(TA)], square interface [CBN(S)], 
and circular interface [CBN(C)]. The size of the interface 
is characterized by the number of carbon atoms at the 



interface, N^dgc, which varies from to 60 in this work. 
Two same phenomena are observed from the two differ- 
ent investigation methods. Firstly, the heat current is 
reduced linearly with the increase of A^cdgc in all sam- 
ples, and the reduction is more than 50% for A^odgc ~ 50. 
Secondly, with same A^odgc, the zigzag interface causes 
stronger reduction of heat current than the armchair in- 
terface. These effects are interpreted from the lattice 
dynamics analysis and the NEGF transmission function 
which are nicely consistent with each other. The lattice 
dynamics analysis finds that the heat reduction is in- 
duced by the localized phonon modes at the CBN inter- 
face. The density of states (DOS) for the localized modes 
increases with increasing TVcdgc and it is larger in zigzag 
interface than the armchair interface, which explain the 
two observed phenomena. Particularly, strong localized 
modes are found in the frequency regions around 300, 
800, and 1400 cm'^ In the other explanation, the NEGF 
transmission function is depressed more seriously with in- 
creasing A^cdgc, and stronger depression is found in the 
zigzag interface than the armchair interface. It also shows 
that the transmission function of phonon modes around 
300, 800, and 1400 cm~^ undergo strong depression for 
all CBN interfaces, which coincides with the former lat- 
tice dynamics analysis. 



II. INTERATOMIC POTENTIAL 

There are three different elements, C, B, and N, in 
present work. We have recently developed a set of Ter- 
soff-f-UOOP potential to simultaneously describe both 
the in-plane and out-of-plane vibrations in hexagonal 
boron nitride sheet^lS This potential is of high efficiency, 
stability, and is particularly suitable for the investiga- 
tion of heat transport as it properly describes the phonon 
spectrum. The interaction between B and N in the CBN 
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FIG. 1: (Color online) Phonon spectrum of graphene cal- 
culated from Tersoff potential (dotted blue line) and Ter- 
soff+UOOP potential (solid red line). Experimental data are 
read from figure 2 of Ref. [T^ 



samples is governed by this set of Tersoff+UOOP po- 
tential. For the interaction between carbon atoms in 
graphene, we can develop a parallel set of Tersoff -1-UOOP 
potential as the graphene and BN sheet have quite similar 
honeycomb structure. The original Tersoff parameters 
for carbon system were developed for sp3 structureSfii so 
they can not be directly applied to sp2-bonding graphene. 
Explicitly, the frequencies of in-plane vibrations from this 
potential are about 60% higher than experimental values. 
To eliminate this flaw, we scale the A and B parameters 
in Tersoff potential by a factor of 0.4, which leads to rea- 
sonable frequencies for in-plane vibrations as shown in 
Fig. [H by dotted line (blue online). However, the result 
also shows that the frequencies of out-of-plane vibrations 
(ZA and ZO) are obviously smaller than the experimental 
value. This disagreement roots in the form of the Tersoff 
potential, where the in-plane and out-of-plane vibrations 
are coupled together. To modify the out-of-plane vibra- 
tions without affecting the in-plane vibrations, we intro- 
duce the universal out-of-plane (UOOP) potential^^ to 
describe the out-of-plane vibrations. The UOOP poten- 
tial has the form V = Co(Ci + C2 cos(0) + C3 cos(2</))), 
where (j) is the dihedral angle. The four optimized pa- 
rameters Ci are -0.588047 eV, 1.000000, 0.814936, and 
0.259776, respectively. This set of Tersoff-^UOOP po- 
tential results in good phonon spectrum of graphene as 
shown in Fig. [T] by solid red line, so it is suitable for the 
study of heat transport in graphene. Note, however, that 
the calculated frequencies of the high energy in-plane op- 
tical modes (LO and TO) do not agree with experiment 
in detail, because these two optical modes are more sen- 
sitive to long-range interaction which is ignored in the 
bond order Tersoff potential. The interactions between 
different elements are also described by the Tersoff po- 
tential with the parameters determined by the combining 
law of multicomponent systems The single parameter 
X introduced by the combining law is set to be 1.0, which 




FIG. 2: (Color online) The heat current v.s A^cdgc for different 
CBN interfaces: (a) and (b) from MD simulation; (c) and (d) 
from NEGF calculation. 

may lead to possible error. However, the error is quite 
small, because x acts as a fine tuning parameter. Fur- 
thermore, the success of the combining law mainly re- 
lies on the Tersoff potentials between individual element, 
which have been carefully fitted to experimental data in 
the above. 



III. RESULTS AND DISCUSSION 

We run MD simulation to calculate the heat current 
J across the CBN compounds. J is calculate from 
J = (Jl — Jfl)/2, where Jl/r are the heat energy pumped 
into system from left /right heat baths. The tempera- 
tures of two baths are Tj^/r = (1 ± a)T with a = 0.1. 
The constant temperature of heat bath is realized by the 
Nose Hoover thermostatj ^^i-^^ A time step of 1.0 fs is used. 
The simulation is run at least for 2 x 10® MD steps and 
will be sufficiently extended to ensure that the steady 
state is achieved. The system is in steady state when the 
difference of heat current from left and right baths is less 
than 2% of the current across the system. We study the 
heat current in the graphene nanoribbon with armchair 
edge. The length of the nanoribbon is 102.24 A and the 
width is 39.36 A. Some C atoms in the center region are 
substituted by B and N atoms. The BN regions in the 
graphene nanoribbon can be of different shapes and sizes. 
Periodic boundary condition is applied in the lateral di- 
rection of the nanoribbon, so that we can avoid the edge 
effect and focus on the effect of the CBN interface. We 
have also performed calculation for graphene nanorib- 
bon with zigzag edge and the same size, where similar 
phenomenon was observed. Fig. [5] (a) and (b) show the 
heat current J across four CBN interface of various sizes. 
With the increase of Ncdgc, a linear decrease of J is ob- 
served for all CBN interfaces. The MD simulation results 
are displayed by points, and the lines are the linear fitting 
of these MD results. For -/Vcdgo ~ 50, the heat current J 
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FIG. 3: (Color online) Vibrational morphology of some se- 
lected localized phonon modes in different CBN compounds 
with Ncdgc ~ 45. The CBN compounds from (a)-(d) are 
CBN(TZ), CBN(TA), CBN(S), and CBN(C). 



is reduced by more than 50%. With same value of A^odgc, 
the CBN(TZ) interface causes the strongest reduction to 
heat current, while the CBN(TA) interface introduces the 
weakest reduction. The heat currents across the CBN(S) 
and CBN(C) interfaces are in the middle of these two 
trigonal samples, where the reduction of J by CBN(C) 
interface is slightly stronger than CBN(S) interface. All 
these phenomena are the same at both 300 K and 1000 K, 
except different value for J. The difference between heat 
current in zigzag and armchair types of CBN interfaces 
increases with increasing A^odgc- A quite large difference 
is desirable in experimental CBN compounds as iVedge is 
usually very large. ^ We note that the quantity Nedge 
is an appropriate choice to distinguish the clus- 
ter types. If the thermal current is plotted as a 
function of cluster area, it is almost a single val- 
ued function of the cluster area, which implies 
that the cluster area also plays a role besides the 
number of edge atoms. 

A straight forward physical interpretation for the heat 
reduction can be made from the lattice dynamics anal- 
ysis. There are some localized modes around the CBN 
interface. These modes do not move and make few contri- 
bution to heat transport. Fig. |3] show selected localized 
modes for all studied CBN interfaces. In these modes. 




FIG. 4: (Color online) Temperature distribution for differ- 
ent CBN compounds: (a), pure graphene; (b). CBN(TZ); 
(c). CBN(TA); (d). CBN(S); (e). CBN(C). Temperature 
decreases from 330 K to 270 K from blue to red. 



only boron and nitride atoms vibrate while all carbon 
atoms do not take part in the vibration. As a result of 
these localized modes, there will be temperature jumps 
at the CBN interfaces as shown in Fig. For different 
types of CBN interfaces with the same iVcdgc, the number 
of atoms per length at the interface is the largest (small- 
est) in zigzag (armchair) type of CBN interface. As a 
result, the DOS of the localized modes in the CBN(TZ) 
is larger than that of the CBN(TA), so the CBN(TZ) in- 
terface induces stronger heat reduction than CBN(TA). 
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FIG. 5: (Color online) Localization properties for phonon 
modes in different CBN interfaces, (a)-(c) are the total DOS, 
projected DOS, and IPR of CBN(TZ), CBN(TA), and pure 
graphene. (d)-(f) are for CBN(TZ) with Arcdge=45, 27, and 
pure graphene. 



The CBN(S) and CBN(C) interfaces are a mixture of 
zigzag and armchair types; thus their reduction capabil- 
ity is in the middle between the two trigonal samples of 
pure zigzag or armchair interfaces. In Fig. [5] (a), (b), 
and (c), we numerically compare the localized modes in 
CBN(TZ), CBN(TA), and pure samples. Panel (a) shows 
that all three samples have similar total DOS. Panel (b) 
is the DOS projected (PDOS) onto the CBN interface. 
There is no CBN interface in pure sample, so PDOS is 
zero in whole frequency range. The PDOS in CBN(TZ) is 
larger than the CBN(TA) due to larger number of atoms 
per length at the CBN(TZ) interface, which confirms our 
above qualitative argument. Panel (c) shows the inverse 
participation ratio (IPR) of each mode. The IPR value 
of the localized modes is about one or two orders larger 
than transmissible modesJ^ There is no localized mode 
in pure sample, and the IPR in whole frequency range is 
very small. In the two trigonal samples, there are lots of 
localized modes. Especially, those localized modes in the 
frequency range of 300, 800, and 1400 cm~^ have very 
high IPR. Panel (b) and (c) together manifest that the 
zigzag interface can introduce stronger heat reduction. 
For a particular sample, the length of CBN interface in- 
creases with increasing -/Vodgo- This leads to more local- 
ized modes at the CBN interface, resulting in stronger 
reduction. Fig. [5] (d), (e), and (f) compare the DOS, 
PDOS and IPR in CBN(TZ) interface with TVedgo = 27, 
45, and pure samples. Panel (e) shows that the PDOS in 
the sample with Nodgo = 45 is obviously larger than that 
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FIG. 6: (Color online) The NEGF transmission function for 
CBN interface of different shapes and sizes. 



of A'odgc = 27, resulting in much stronger reduction to J. 

In the above, we have run MD to simulate the heat 
current adjusted by various CBN interfaces, and pro- 
vided an explanation by means of the lattice dynamics 
analysis. However, there is no phonon-phonon scatter- 
ing effect in the localized modes, while this effect is con- 
sidered in MD simulation. To reinforce our conclusions 
from MD simulation and the explanations from localized 
modes, we repeat the whole calculation using the ballistic 
NEGF approach. 18 The ballistic NEGF do not take into 
account the phonon-phonon scattering. But it is an exact 
treatment for the scattering at the CBN interface, which 
actually reflects the localization property of the localized 
modes. If the NEGF leads to the same conclusions as MD 
simulation, then we can safely announce that the phonon- 
phonon scattering in MD simulation has no effect on heat 
reduction. Fig. [2] (c) and (d) show the heat current cal- 
culated from NEGF, which are consistent with MD sim- 
ulation results. Firstly, J decreases almost linearly with 
increasing A^cdgc for a particular sample. For A'edgc ~ 50, 
the heat current is reduced by more than 50%. Secondly, 
for the same N^dgc, the CBN(TZ)/CBN(TA) interfaces 
cause the strongest /weakest reduction of heat current. 
Heat reduction from the other two interfaces are be- 
tween these two trigonal interfaces. These two results 
are the same at both 300 K and 1000 K. It should be 
note that the value of J from NEGF approach is almost 
twice of that from MD simulation at 1000 K, because 
of strong phonon-phonon scattering in MD simulation 
at high temperatures. The value of J from two differ- 
ent methods is close to each other at 300 K, where the 
phonon-phonon scattering is weak. We have obtained the 
same conclusions from MD simulation and NEGF calcu- 
lation, although the values of J are quite different from 
these two approaches. We further show the transmis- 
sion function in the whole frequency range from NEGF 
method in Fig.[B] It shows that the transmission function 
in CBN(TZ) is smaller than that of the CBN(TA) with 
same A^cdgc- The transmission function is smaller in sam- 



5 



I 10'' r 












zigzag 






armchair _ 
















1000 2000 





%Q0 600 900 1260 1'sod 1800 
T(K) 



FIG. 7: (Color online) The interface thermal resistance v.s 
temperature for zigzag (solid line, green online) and armchair 
(dotted line, red online) CBN interfaces. Inset shows the in- 
terface thermal resistance (in logiox scale) in whole temper- 
ature range. 

pies with larger iVodgc, indicating stronger heat reduction 
in these samples. A general feature is that the transmis- 
sion functions for phonon modes with frequencies around 
300, 800, and 1400 cm~^ are depressed seriously by CBN 
interfaces. This feature agrees quite well with the strong 
localized modes in these frequency regions as we have 
discussed previously. These results further confirm that 
the reduction of heat current is because of the localized 
modes at the CBN interface. Fig. [6] also shows that high 
frequency modes, which have increased contribution to 
heat current at higher temperatures, are suppressed more 
effectively. As a result, the reduction is much more effec- 
tive at high temperatures in Fig. [2l as indicated by the 
larger slope of the fitted line at 1000 K than that of the 
300 K. 

The NEGF approach can provide us another more im- 
mediate evidence-interface thermal resistance (ITR)-for 
the different heat reduction in zigzag and armchair type 
CBN interfaces. We apply NEGF to study heat transport 
across junctions formed by graphene on left part and BN 
sheet on the right part, with the CBN interface at the 
middle. NEGF calculates the thermal conductance (a) 
of the junction. The thermal resistance can be obtained 
through R — s/a, where s is the cross section of the junc- 
tion considering 3.35 A as the thickness. The ITR can be 
estimated by ITR= R — Rq, with Rq as the thermal resis- 
tance in the pure graphene sample. Fig. [7] compares the 
ITR for zigzag and armchair CBN interfaces. The ITR 



in zigzag CBN interface is about 7 x 10~-^°m^K/W at 
room temperature and decreases to 3.5 x 10^^°m^K/W 
at higher temperatures. The value of ITR in armchair 
CBN interface is only about half of the zigzag CBN in- 
terface in whole temperature range. As a direct result 
from its higher ITR, the zigzag CBN interface can in- 
troduce stronger heat reduction than the armchair CBN 
interface. 

IV. CONCLUSION 

To conclude, this work studies how the heat current 
is manipulated by the CBN interface of various shapes 
and sizes. The interatomic interaction is described by a 
set of Tersoff-|-UOOP potential, which is efficient, sta- 
ble, and suitable for the study of heat transport in the 
sp2-bonding system. The calculation is accomplished 
by using two distinct approaches: MD simulation and 
NEGF calculation. Both methods find same phenom- 
ena that the heat current is reduced linearly with the 
increase of A^odgo, and the reduction can be more than 
50% for iVedgo ~ 50. It is also shown that the zigzag 
CBN interface introduces stronger reduction to heat cur- 
rent than the armchair CBN interface. These observa- 
tions are explained from two different points of view, 
which are eventually consistent with each other. Both 
explanations interpret that the localized phonon modes 
at the CBN interface play a key role in the heat reduc- 
tion and the phonon-phonon scattering has no effect. A 
more immediate evidence from NEGF for the stronger 
heat reduction in zigzag CBN interface is its stronger 
ITR, which is nearly double of that in the armchair CBN 
interface. 

The graphene and BN sheet are two promising nano- 
materials with plenty of advance properties. They have 
similar honeycomb structures with close configuration 
parameters fi2i^ Existing works have demonstrated that 
the compounds of these two materials can upgrade the 
electronic properties in graphene and particularly tun- 
able electronic band gap can be realized. The practi- 
cal significance of present work is to provide a guideline 
for the experimentalists to realize tunable heat transport 
capabilities through shaping the CBN interfaces during 
the synthesis of CBN compounds. So that higher heat 
transport capability can be achieved for rapid electronic 
devices and lower heat transport capability is obtained 
for thermoelectric nano-devices. 
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